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APPROXIMATION BY PLANAR ELASTIC CURVES 


DAVID BRANDER, JENS GRAVESEN, AND TORE BJERGE N0RBJERG 


Abstract. We give an algorithm for approximating a given plane curve segment by a pla¬ 
nar elastic curve. The method depends on an analytic representation of the space of elastic 
curve segments, together with a geometric method for obtaining a good initial guess for the 
approximating curve. A gradient-driven optimization is then used to find the approximating 
elastic curve. 


1. Introduction 

An Euler elastica or elastic curve is the solution to the variational problem of minimizing 
the bending energy, the integral of the eurvature squared f K:(s)^ds, among eurves of a given 
length with fixed endpoints and with the tangents preseribed at the endpoints. All solutions 
to this problem were deseribed by Euler [4] in 1744, and the eurves ean be parameterized in 
terms of elliptie funetions. 

The energy minimizing property qualifies fhe elasfiea as a mafhemafieal model for fhe 
shape assumed by a fhin inexfensible rod when eonsfrainfs are plaee only af fhe endpoinfs. 
This shape appears nafurally in eerfain manufaefuring seenarios: for example in a eonsfrue- 
fion made from fhin, flexible sfrips of wood or similar, fhe shape of eaeh sfrip, befween fhe 
fixed poinfs, is an elasfiea. In anofher applieafion, a fhin mefal blade ean be healed and used 
lo euf polyslyrene for arehileelural formwork, fhe blade ehanging ils shape during fhe mo- 
lion. This permils fhe eonslruelion of eurved geomelries polenlially mueh more eheaply lhan 
wilh fhe alfernafive of numerieally eonfrolled milling; Ihis arliele is pari of a larger projeel 
dediealed lo fhe developmenf of Ibis so-ealled robotic hot-blade cutting leehnology [3]. Here, 
loo, fhe shape of fhe blade is an elasfiea. 

The presenl slandard represenlalion for eurves and surfaees in eompuler-aided design 
(CAD) sysfems is rational splines, lhal is, pieeewise ralional funelions. In arehileelure Ihere 
is lypieally a seeond slep after fhe inilial eoneeplual design, whereby fhe CAD model is 
adapled slighlly for praelieal realization wilh respeef lo Ihe ehosen manufaefuring proeess, 
and Ihis is ealled rationalization. The problem of rationalizing a slandard CAD design for a 
eonslruelion melhod involving elaslie eurves enlails Ihe approximation of a ralional or poly¬ 
nomial spline eurve by an elasliea, and il is Ibis problem lhal motivates us here. Il is worlh 
noting lhal Ihe term “spline”, in pre-CAD years, referred lo a Ihin slrip of wood used lo draw 
eurves lhal interpolate a given sel of poinls smoolhly. Belween Ihe interpolation poinls, Ihese 
slrips assumed Ihe shapes of elasliea, or ralher planar elasliea, beeause Ihe slrips were laid 
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on a flat surface. Thus, for the type of rationalization referred to above, one could say that 
the task is, somewhat ironically, to approximate a digital spline with an analogue spline. 

With this in mind, we consider here the problem of approximating an arbitrary curve by a 
planar elastica. 

Precursors: Quite distinctly from the applications mentioned above, elastic curves are a 
good candidate for the choice of functions used in geometric modeling because their curva¬ 
ture minimizing property makes them optimally faired. The idea of using them as mathe¬ 
matical splines - here meaning piecewise smooth functions - has been discussed by many 
authors: for example the references [1], [2], [5], [6], [11], [12] and [13] constitute a repre¬ 
sentative, though not an exhaustive, list. These works are concerned either with the problem 
of how to compute the elastic curve satisfying various constraints respecting placement of 
endpoints, end-tangents, lengths etc., or with such problems as the existence of an elastic 
curve interpolating a given set of points. As a testament to the importance of these curves, 
one finds that many aspects of the theory are re-derived independently several times. 

For our purpose, however, the main point is that previous work on elastica and their appli¬ 
cations in spline theory does not consider the problem of approximating an arbitrary curve 
by an elastica. 

Outline of this article: We first make use of the analytic solutions, which were derived in 
1880 by L. Saalschiitz [15], to give a representation of the space of elastic curve segments. 
The space depends real analytically on seven control parameters, and one can therefore, in 
principle, use gradient driven optimization software to achieve an approximation. The chal¬ 
lenge is that the result of this non-convex nonlinear optimization depends very much on the 
initial guess, as illustrated in Figure 1. 






Figure 1 . The blue curve is to be approximated by an elastica segment. 

The green curves show the result of IPOPT optimization with different initial 
guesses. In the third case the optimization terminated before an extremum 
was reached. 

The optimization approach can only deliver a useable approximation algorithm if a means 
of choosing a good initial elastic curve is found. Our solution is to use the fact that the 
curvature function for an arc-length parameterized elastica is affine in a certain direction. We 
first show, in Section 3, how to use this fact to recover the seven control parameters - in a 
numerically stable manner - from a given elastic curve segment. We can then (Section 4) 
apply essentially the same procedure to an arbitrary plane curve to obtain a canonical 
elastic curve segment that, with respect to the global geometric characteristics selected for in 
the procedure, is close to the given curve. 
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The canonical initial guess can then he taken as input for an optimization package to 
compute a good approximating elastic curve for a given curve, provided a good solution 
exists. We illustrate the results of the overall algorithm on a sample of Bezier curves - some 
of which are close to elastica and some of which are not - with endpoints free (Figure 7) and 
fixed (Figure 8 ). 

Finally, in Section 5, we discuss applications of this method to the problem of approximat¬ 
ing curves hy piecewise elastic spline curves and ongoing work on applications in manufac¬ 
turing. 


2. Planar Elastica 


2.1. Euler-Lagrange equation. Here is a brief sketch of the equations defining planar elas¬ 
tica. More details, background and references can be found in [ 8 ]. Let 7 : [0,£] —)• be 

a plane curve segment parameterized by arclength. Let G (5) denote the tangent angle, de¬ 
fined by the equation 7 ( 5 ) = (cos 0 ( 5 ), sin 0 ( 5 )). Then a curve segment of length £ starting at 
(xo,yo) and ending at satisfies = xq -|- /gCos Gds and = yo + fo sin 0 d 5 . 

Lef K denote the curvature 6(s). An elastica is a minimizer, among curves with the same 
endpoints and end tangents, of the bending energy ^ Jq fc( 5 )^d 3 '. Suppose 7 is an elastica 
from (xo,yo) to (x£,y^) with angle function 0 , and consider the perturbed curve 7 , with angle 
function 6t(s) = 6(s) + ty/(s), where i/r is a differentiable function with y/(0) = y/(£) = 0. 
Applying the method of Lagrange multipliers to the bending energy, we set: 


<^( 7 ) = 2 J di-hAi^xo-h J cosGds-Xi^+h(yo + sin0d5-y^^ 

and we require that 

d^(7j I 


0 = 


dt 


f =0 


fi /d^0 

~~ “TT + ■^1 0 — ^2 cos 0 ) di. 


/o 




Since i/r was arbitrary, it follows that 0 satisfies the Euler-Lagrange equation 


( 2 . 1 ) 


—^ -|- Ai sin 0 — A 2 cos 0 = 0 . 
ds^ 


Setting (Ai, A 2 ) = A (cos 0, sin ^), with A > 0, this becomes 0-|-Asin(0 — = 0. Note that 

A = 0 if and only if K is constant, i.e. the curve 7 is either a straight line segment or a piece of 
a circle. If A / 0, set 7(5) = s/XR-^jis / \/A), where is the rotation by angle 0. Then yis 
also an elastica with tangent angle 6(s) = 6 (s/VX) — 0 satisfying the normalized pendulum 
equation 6" = — sin0. In summary: 


Theorem 2.1, Up to a scaling and rotation of the ambient space, all arclength parameterized 
elastica 7 : [ 0 , 1 ] —)• with non-constant curvature K, can be expressed as: 7(5) = 7 ( 0 ) -|- 

/o (cos 0 (t), sin 0 (t))dt, where 


(2.2) 


0 = — sin0. 
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2.2. Parameterizations of the space of elastica. We now want to define some suitable con¬ 
trol parameters to describe an arbitrary elastic curve segment. Essentially, the parameters 
need to specify which solution to (2.2) is involved, the start and endpoints on the solution 
curve in question, and a rotation and scaling. 

The elastic curves can be expressed in closed form via the elliptic functions (see Appen¬ 
dix A). The formulas can be found in Love [10]. There are two classes of elastica: curves 
with inflection points (i.e. points where 6=0) and curves without inflections. 

Basic elastica. The solution to (2.2) starting at (0,0) with initial angle 6 (0) = 0 and 6 (0) > 0 
is 

= {2E{s,k)-s, 2k{l-cn{s,k))), 

where k = 6{0) /2. For ^ G [0,1) we get inflectional elastica; for ^ > 1, we use the extended 
elliptic functions defined in Appendix A to obtain elastic curves without inflections. We 
reserve the name for these basic elastica. Figure 2 shows elastic curves for different 
values of k. All elastica are obtained by scaling and rotating these curves. The periodicity of 
the curves is given by: 


i;,{s + AK) = + {2E{AK)-AK, 0), 


where K the quarter-period defined in Appendix A. 


k=Q 


k=0.3 





6mm 




k = 10 



fc=0.937 


fc—>CXD 



Figure 2. Elastica for different values of k. The curves are scaled to a 
uniform “height”. 
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General elastica. To parameterize an arbitrary segment of an elastiea, we ean ehoose a seg¬ 
ment of a basie elastiea by ehoosing k, a starting point and an endpoint sq + £, where 
f'GM\{0}. It will be eonvenient to remove the dependenee on sq and £ from the domain 
and inelude them in the parameterization. We parameterize the elastiea segments on the unit 
interval, setting s = so + £t, with t G [0,1]. The new eurve parameter t is not unit speed. Fi¬ 
nally, any elastiea segment ean be obtained by introdueing a sealing faetor 5 > 0, a rotation 
by an angle (j) G (—71, Tl] and translation by a veetor (xo,yo)- We thus have a standard elastie 
segment parameterization 

7(Uo,^,5,(^Ao,yo)(^) = +^0 + (-^0, f'o), I £ [0)1]- 

It depends on seven eontrol parameters, but we will usually omit the subseript. Sueh a eurve 
has eonstant speed \l\S and length L = 

Remark. If £ < 0, the orientation of the eurve is ehanged. In the infleetional ease, the elastiea 
with opposite orientation ean be obtained by a rotation by n. In this ease we may therefore 
assume £>0 without loss of generality. For elastiea without infleetions, however, a segment 
with £ <0 eannot be deseribed as a segment with £>0. One ean instead reverse the direetion 
of the parameterization for that ease, and so all eases ean be handled with the assumption 
£> 0 . 

For any elastie eurve 7 of the above type, the eurve KO = 7(s) Letting d 

and d denote the angle funetions of and 7, respeetively, we have d(t) = 0(^0 + y) + 0? 
and thus 

S"(t) = ^0(50 + 1) = -^sin0 (^o + y) = -^sin(0 (t)- 0 ), 
so the angle funetion for 7 satisfies (2.1) with 

(2.3) (Ai, A2) = ^(eos0, sin0). 

We will also use the faet that the eurvature for the elastiea 7^^ soiS(t>xo yo) 

(2.4) K(t) = Ycn(so + £t). 

3. Finding the control parameters of an elastic curve segment 

We first deseribe a way to ealeulate numerieally the eontrol parameters of a given planar 
elastie eurve segment. In the next seetion the same reeipe will be applied to an arbitrary 
planar eurve to obtain a eanonieal first guess for an approximating elastie eurve. The main 
idea is to exploit the faet that the eurvature of an elastiea is an affine funefion of the distanee 
along a speeial direetion. 

Let X : [a, fi] —)• be an elastie eurve parameterized by arelength. As for any planar eurve, 
we ean write the tangent and the normal as t = (eos0,sin0), n = (—sin0,eos0), and we 
have the Frenet-Serret equations 

dt d0 d^t (fe 2 

The tangent angle 0 must satisfy the Euler-Lagrange equation (2.1) for some Lagrangian 
multipliers , X 2 to be found. 
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Let u denote the projeetion of x onto the line spanned by (A2, —Ai), i.e., 

1 . Xix-X\y 

U = ^(A 2 , -M) ■ {x,y) = --, 

where A = || (Ai, A2) || = Setting ^ = 0 in (2.3), we find that the veetor (A2, —Ai) points in 
the downward direetion in Figure 2. It follows that u is bounded and periodie in s. Moreover, 
we ean write the Euler-Lagrange equation as 0 = A d, so we have 


(3.1) 


K=-^=Xu + a = X2X — Aiv- 
di 


a, 


whieh is to say that the eurvature is an affine funefion of u. 

In order fo find Ai, A2 and a in a numerieally sfable manner, we solve fhe above equation 
in fhe leasf squares sense, i.e., we eonsider fhe quadrafie minimization problem 

rb 


f 2 

limize / (k: +Aiy — A2X — a) d^, 
,X2,ot Ja 


minimize 
A, 


whieh leads fo fhe following linear sysfem 

t^y^ds -l^xyds -l^yds 


(3.2) 


Ja ' 

-f^xyds J^x^ds 




J^XKds 


Eel 6u denote fhe angle befween fhe fangenf veefor t and fhe n-axis (see Eigure 3). We 
have 

dx du 
(E “ (E ’ 

sin0„ = y(Ai, A2) -t, 


(3.3) 


eose„ = ^(A2, -Ai)-t= ^(A2 , -Ai) 


and henee 

or equivalenfly 
(3.4) 


dsine„ 

du 


A 

di dsin0„ 
du di 


1 


a ddu , , 

eos 6u — = ic = A u + cc, 
cos du di 


1 


P{u) := sin0„ = -Au^ + au + fd . 


We solve fhis equafion wifh respeef fo j3 in fhe leasf squares sense and obfain 

1 


(3.5) 


j8 = - / ( sin 0„ — ^ A — a M ) di, 
L Ja \ 2 


where L = b — ais fhe lengfh of fhe eurve x. 

Eor an elasfiea x(5) = SRipl^ki^/S) + {xo,yo) we have (Ai,A 2 ) =5'^^(eos(/),sin(/)) and k{s) = 
{2k/S)cn{s/S). Subsfifufing fhese info fhe definitions u = (A2X —Aiy)/A, a = K — Au and 
sin0„ = (l/A)(Ai,A2) -t we have 

u = —25^(1—en(5/5'))+xo sin 0—yo cos 0, 
sin6„ = 2dn^(5/5') — 1, 

a = 2kS^^ cn{s/S) — Au = 2k/S+{yocos(j)—xosin(l))/S^ 


(3.6) 
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Figure 3. A segment of an elastica in the plane, showing 0, 0, the M-axis 
and the angle 


Then the equation j8 = sin 0„ — Am^/2 — au becomes: 


(3.7) 


^ = 1 + 


xq sin (j) —yo cos (j) 
2S2 


(aq sin ^ — yo cos 0 


It follows from (3.4) that all points on the elastica correspond 
of the polynomial P is between — 1 and 1, and hence 

-a-5 -a+ 5- 


u G 


A 


-AkS). 

to M-values where the value 


where 5 = \J —2A(j8 — 1). 

If the elastica has an inflection, there must be some «*, such that K = Xu^, + oc = 0, but this 
means that u* is a minimizer for P{u), and thus its minimum must lie in [—1,1]. Moreover, 
the inflectional elastica has points where sin 0 m = 1 (which happens twice per period), but no 
points where sin 0„ = — 1 (see Figure 4). Hence u runs through all of the interval where P{u) 
is less that 1; in other words, Umin and Umax are exactly the endpoints of the above interval. 

For the elastica without inflections, the tangent makes full rotations, so sin 6u takes both of 
the values ± 1. Hence, in this case, we must have 


.^min) ^max 


-a-5 -a-5+ 

A ’ A 


or 


[Wminj ^max 


—(X 3 

A ’ A 


where 5+ = — 2A(j3 + 1); these are the two cases corresponding to £ < 0 and £ > 0, re¬ 

spectively. We can thus determine whether the elastica has inflection points based on whether 
the minimum for the polynomial P{u) is smaller or greater than —1, see Figure 5. In fact. 



















DAVID BRANDER, JENS GRAVESEN, AND TORE BJERGE N0RBJERG 




Figure 4. The inflectional elastica have points where sin0„ = 1, but not 
— 1. For the non-inflectional elastica sin0„ takes both the values ±1. 


from (3.6) and (3.7) we have 

or equivalently 
(3.8) 


^ 4^2 

a2-2A(^-l) = ^ 


k = 


y/a^-2X{p-l) 

iVI 


so we can find S, (j) and k from Xi, X 2 , oc and jS. 


Remark. The above formula also holds if ^ < 0 and so does the expression for j3 in control 
parameters. The expressions for fc, X\, X 2 and a simply change sign in this case. 




Figure 5. The parabola (3.4). To the left in the case of an elastica with 
inflection points, to the right without. The blue and red part corresponds to 
points with negative and positive curvature, respectively. 

We still need to recover sq and i. We have 

u = —2kS (1 — cn (^0 + I)) +xosin0 —yocos(j), 





















APPROXIMATION BY PLANAR ELASTIC CURVES 


9 


and since 


we get 


Mmax = =xosin{^)-yocos{^), 


A(u) = Mmax — u = 2kS { \ — cn^) , 


(3.9) cn(5,^) = l-^. 

If we consider the unbounded complete elastica, then u oscillates between Wmin and Umax and 
we can divide the elastica into segments where u is monotone, each with length equal to a 
half period 2KS. 

We first consider the case of an elastica with inflection points (i.e. k < 1). Here we 
have cn(5,fc) = cos (am(5,^)). If the start point xq = x(a) is on segment number 1 and u is 
decreasing here, then 

am{sQ,k) = arccos (1 - , 


if n is odd, 
if n is even. 


and if the end point xi = x{b) is on segment number n, then 


am(5i,^) = 


[n — 1) 71 +arccos |^1 — 
nn — arccos 1 —++ 


2ks ; ’ 


2kS 1 ’ 



In all cases we have 

Si = F{sLm{si,k),k), f = 0,1, 

and l = s\ — So- 

In the case of an elastica without inflections points (i.e. ^ > 1) we need a little work to find 
am. We have 

sn{s,k) = ^ ~ k 

and 

sn{s,k) = 

If u is decreasing on segment 1 then 
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and if we have n segments 


am ( ^, ^ ) = < 
k 


j" 71 + aresin ^^k— , if n is odd 


171 — aresin^ 


if n is even. 


If u is inereasing on segment 1 then 


am {— ,k] = 71 — aresin< 
k 


Ia{uq) ( A{uq) \ 
AS ) 


and if we have n segments 




a.m{j,k] = { 


j" ^71 — aresin ^^(|i) , if n is odd, 

jTT + aresin ^[k - , 


if n is even. 


Finally, we find the 5-values using the ineomplete elliptie integral 

Si = {F (am [ksi, i = 0,1, 


and (. = s\ —50. 

Remark. If we have a negatively eurved noninfleetional elastiea (i.e. I < 0), we ean reverse 
the parameterization, find fhe elasfiea, and inferehange (50,5i). 


We now have a sealed and rofafed elasfiea segmenf, 7o = /(jt ^ 5 ^ 0 o)’ 
is fo find fhe final franslafion (xo,yo)- This is done by solving fhe equation 

= 7o('^) + (-^o,7o), 

in fhe leasf squares sense. The solution is 

(3.10) (xo,yo) = (x(5)-7o(5))d5. 


4. Approximating a plane curve by a planar elastica 


We are now given a eurve x: [0,1] —)■ M^, nol neeessarily elasfie and nol neeessarily pa- 
rameferized by arelengfh. The arelengfh is given by 

5(0= riix'(T)iidT, 

Jo 

and fhe lengfh of fhe eurve is L = 5(1). We wanf fo approximate fhis eurve by a pieee of an 
elasfiea. We do fhis by minimizing a suifable disfanee, sueh as fhe L^, or disfanee, 
over fhe eonfrol parameters p = {k,so,(,S,(j),xo,yo)- In the ease of L? fhe problem is 

minimize ,^(fc,5o,f?,5',0,xo,yo), 

k,so,e,S,(j>,xo,yo 


where 



x'(f)||dt. 
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If we want the elastie eurve to satisfy further eonditions, such as having the same endpoints 
and/or end tangents as the original curve, we can include these in the optimization problem. 

For the optimization, we have used the gradient driven tool IPOPT [17], so we need the 
first and second order partial derivatives of with respect to the control parameters, which 
are 


dpi 

dpidpj 



^ f 

^P‘ \ L ) 


x'(t)||dt, 





^Pj \ L ) 



dpidpj 



x'(t)||dt. 


See Appendix B for a list of specific derivatives. 

The optimization problem is non convex and the result depends on the initial guess (see 
Figure 1). A canonical geometrically plausible guess is obtained from a generalization of the 
procedure of Section 3 to the case of an arbitrary input curve, which we will now describe. 

We find Xi,X 2 , oc as before, by solving (3.2), noting that = Jq f{t) ^dt. This 

gives us the scaling and rotation of the elastica. We can judge the success by calculating the 
normalized residual 


Ri = 




— dt. 
dt 


Similarly j3 can be found by (3.5), where sin0„ is given by (3.3), and we can calculate the 
normalized residual 



Remark. Another possibility is to forget that we know A and a and solve (3.4) with respect 
to A, a, and jS, but in the few cases we tried this, the results got worse. 

We know that sin takes values in [—1,1], and since fd is chosen to minimize the distance 
between sin and the polynomial P{u) = | A + a « + j8, the latter must be less than 1 for 
some M-values, so the number 5 = — 2A(j3 — 1) is well-defined. We can thus determine 

whether the elastica has inflection points and we can determine the parameter k from (3.8). 

At this point we need to take into account the fact that the input curve is not necessarily 
an elastica. For an elastica, we could easily count the oscillations, but for an arbitrary curve 
there may be oscillations of different sizes. We find the curve segments where u is monotone, 
but we only count such a segment as an oscillation if it has some minimal height: we have 
used half of the difference Umax — Umin (as defined in Section 3) as this minimum. Moreover, 
the right hand side of (3.9) need not be between —1 and 1, or, in the noninflectional case, 
between 1 — 1 and 1. We have circumvented this problem by replacing too small values 
by —1 (or y^l — 1/k^) and too large values by 1. The two issues are illustrated in Figure 6. 
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segment), it is simply cut off in these regions. The resulting elastica is shown 
in green. On the left all oscillations of the input curve are counted, on the 
right the two very small ones are ignored. 

We can thus find and We can judge the validity by calculating 














Figure 7. Examples of cubic Bezier curves approximated with elastica. 
The solid blue line is the Bezier curve and the dashed red line is the ini¬ 
tial guess for an approximating elastica. The solid green curve is the best 
approximating elastica found with IPOPT optimization. 
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Figure 8 . The solid cyan curve is the best approximating elastica having 
the same endpoints as the original curve (blue). The green curve is the output 
curve with free endpoints, as in Figure 7. 

We finally determine the translation by (3.10) and define the residual as 



where po is the vector of control parameters found by the above procedure. 

We have tested the procedure on a selection of cubic Bezier curves, displayed in Figures 7 
and 8. In Table 1 we have reported the residuals. 

5. Conclusions, discussion and related work 

When combined with a suitable segmentation, we have found that the method described 
here gives an effective algorithm for approximating curves by piecewise elastic curves. The 
degrees of freedom allow piecewise elastica, with continuity at the joins if end-points are 
fixed, and with continuity if end-points are allowed to move. The method is incorporated 
in on-going work on approximating surfaces by segmented surfaces, with the segments swept 
out by elastic curves (Figure 9), for the purpose of manufacturing architectural formwork by 
robotic hot-blade cutting [3], [16]. 

Our choice of parameters allows us to work with analytic expressions of elastic curves. 
The advantage is that when the seven parameters are known any subsequent calculation is 
accurate and easy to perform. The disadvantage is that the geometric meaning is not obvious 
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Table 1. The first column refers to the examples in Figure 7, the next three 
report the residuals Ri, R 2 , and R^, in the approximation process. Next, 
we have the normalized L^-distance, R 4 = for the initial guess 

and then R 4 , the gradient norm || V#'|| and the number of iterations for the 
optimized elastica without and with endpoint constraints. We use the same 
initial guess, whether we constrain the endpoints or not. 



R\ 

Ri 

/?3 

t^4(Po) 

^4(popt) 

l|V^(Popt)|| 

ttiter 

^4(pSpt) 

U iter* 

1 

0.46 

0.14 

0.0 

0.0097 

0.0080 

2.2e-08 

35 

0.0081 

10 

2 

0.45 

0.23 

0.52 

0.048 

0.038 

5.1 

lOOOt 

0.063 

lOOOt 

3 

0.14 

0.020 

0.19 

0.077 

0.0018 

6.0e-09 

26 

0.0023 

16 

4 

0.68 

0.17 

0.14 

0.022 

0.012 

5.9e-09 

8 

0.014 

9 

5 

0.65 

0.27 

0.14 

0.031 

0.018 

5.2e-09 

9 

0.021 

8 

6 

0.099 

0.025 

0.048 

0.0036 

0.0011 

8 .1e-ll 

24 

0.0015 

14 

7 

0.064 

0.0044 

0.048 

0.0010 

0.00032 

1.7e-10 

20 

0.00041 

104 

8 

0.27 

0.069 

0.14 

0.011 

0.0032 

1.5e-09 

199 

0.0046 

too 

9 

0.017 

0.00033 

0.0 

0.00013 

0.0012 

0.053 

1000 

5.1e-05 

15 

10 

0.020 

0.00062 

0.0 

0.00012 

9.9e-05 

1.2e-09 

165 

0.00011 

29 

11 

0.33 

0.014 

0.0 

0.0050 

0.0017 

8.3e-ll 

76 

0.0020 

11 

12 

0.36 

0.10 

0.19 

0.015 

0.0041 

7.3e-09 

178 

0.0053 

83 


flPOPT terminated because iteration count reached maximum (which was set to 1000). 


for all seven parameters. A more geometric set of parameters is the length and the end points 
and tangents; the problem with that alternative is that we would have to solve a nonlinear 
boundary problem to determine the curve and, more severely, the solution is not unique. For 
example, if we imagine that we rotate one of the tangents through 2 % and follow a continuous 
family of elastic curves then the curve at the end will be different from the curve at the start. 

Figure 9. Surface rationalization for hot-blade cutting. Left: a CAD sur¬ 
face patch is approximated by a family of elastic curves. Right: a surface 
segmented into elastica-foliated pieces. 

Note that all the examples in Section 4 have, for convenience, been computed using the 
distance. Given that our initial guess method is based on the curvature, an norm may be a 
more natural choice, but in practice more complicated. 

Another approach to the problem is to minimize / Hy —x|| ds-|-j 8 / K^ds. Here it is not 
clear how to choose j3 . If j3 = 0 then we obtain x = y and in the limit where j3 —)• 00 we obtain 
the best line segment approximation to 7 . In order to make a purely numerical approach work 
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we need to have a good approximation to an elastic curve at all times while we minimize the 
distance to the target curve 7. This means we will have to solve a nonlinear equation at each 
step of the optimization and also find the sensitivities of this solution. 

It should be pointed out that the method for obtaining the initial guess (Section 3) is only 
useful for a curve segment that is not too far from some elastic curve segment. Therefore, the 
practical use of this algorithm requires that a curve first be segmented into suitable pieces. 
There are several ways to approach this: for example (a) apply the initial guess method to the 
curve, (b) measure the distance between the resulting elastic segment and the original curve, 
(c) if the distance is too large, divide the curve into two and repeat. An example with various 
segmentations is given in Figure 10. 



Figure 10. Approximations of a more complex curve (blue) by, in order, 
one, two and four tangent continuous elastica segments. Both the endpoints 
and endtangents of the target curve are matched. 


In physical applications such as hot-blade cutting, it will also be necessary to consider the 
stability of the elastic curve segment obtained from this method; that is, small perturbations 
of the length, endpoints and tangents of the curve should correspond to small changes in the 
solution shape. See, for example, [9] for a study of stability. Constraints such as demanding 
that the curve segments have no inflection points, or an upper bound on the curvature can be 
added to the procedure to ensure stability. 


Appendix A. Elliptic functions 


We list some details of the Jacobi elliptic functions for convenience and to fix conventions. 
Let (0,1). The elliptic functions sn, cn and dn with (elliptic) modulus k are defined as the 
solutions to the system of differential equations: 


sn'(M) =cn(M)dn(M), 
cn'(M) = — sn(M)dn(M), 
dn'(M) = —k^ sn{u) cn{u) 


sn(0) = 0, 
cn(0) = 1, 
dn(0) = 1. 


The complementary modulus k' G [0,1] is defined by k^ + k''^ = \. We have the identities: 


9 7 9 9 9 

(A.l) snM + cnM = l, dnM + ^snM = l, 


dn^ u — k^ cn^ u = k'^. 
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The elliptic functions can be expressed in terms of integrals of trigonometric functions as 
follows. Define the (elliptic) amplitude am as: 

am(t) F{<^)= j ^ -Au. 

•^0 V 1 — sin^ u 

Then 

sn(M) = sin(amM), cn(M) = cos(amM), dn(M) = \j 1 — IF-uv?'(amu). 

Elliptic integrals. The integral F((p,k) given, for each k, by the formula F((p) above, is 
called the incomplete elliptic integral of the first kind. We define the incomplete elliptic 
integral of the second kind by 

E(^,k)= / dn^(M,^)dM. 

Jo 

The complete ellipfic infegrals of fhe firsf and second kind are respecfively fhe funcfions 
F{k)=F{^,k) andE(k) =E(f,k). 


Addition formulas and periodicity. The elliptic functions satisfy the addition formulae: 

snMcnvdnv + snvcnudnM 


sn(M + v) = 
cn(M + v) = 
dn(M + v) = 


l—k^ sn^ u sn^ v 
cnucnv — snusnvdnudnv 
I —k^ sn^ u sn^ v 
dnudnv — snusuvcnucnv 


l—k^ sn^ u sn^ v 


We define the quarter period K by: 

am(.^f) = 


Tl 

2 ’ 


i.e. K = F 


so that sn.^f = 1, cn.^f = 0 and AnK = k'. Then one obtains the periodicity: 

sn(M + 2.k') = — snn, cn(M+ 2.^f) = — cnu, dn(M+ 2.^^) = dnu. 


Extension of k-domain. The elliptic functions, as we have defined them, are only valid for 
k G [0,1]. However, by analytic continuation (see e.g. Lawden [7]), the domain of k may be 
extended. For k> the following identities hold for all n G M: 

sn(M,k) = \ sn(ku, ^), 

cn{u,k) = dn{ku, ^), 

dn{u,k) = cn{ku, 

E{u,k) = kE{ku, ]:) + m(1 —k^). 

Observe that K{k) ^ o° as k ^ 1, so we cannot extend K continuously. We choose the 
extension 

K{k) = T,KCt), k>\, 

which ensures that the period of cn is always 4K. We stress that this is not the analytic 
continuation of K, which in fact takes non-real values for k > 1. 
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Appendix B. Derivatives 

In this section, we list the derivatives of the basic elastica using the shorthand notation 
S = sn{s,k), C = cn{s,k), D = sn(5,^), E = E{s,k). 

We have 




The derivatives with respect to s follow from the definitions of the elliptic functions. 




2D2-1 

IkSD 


ds^ 


^{s,k)=2kC 


-2kSD 


-2kSD\ 
2D2- 1 ) • 


The derivatives with respect to k of sn, cn, dn and E can be found in [7]. From these, with 
repeated use of (A.l), one can find 

— r( ^ 

dk ^ ~ k'^ U'2 + C{k^-D^)-SD{E- sk'^) 


dsdk 


Lf- 

2 


i{s,k) =2 {sD-C(£-sC-)) , 


+ 


2SDC (d2 - k^E^ + k'^ {s^k^ -{E-sY-\)) 
i ((1 - 2k^S'^){E -s){2sk^ +E-s)C + DS{sk'^ -E){Ak^C^ + k'^)) 
{E - s) (C2 + d2 - 4C2d2) + 2sk^{2S^ sk'^ 

-skk'^DS + s^k^C + kCS^ (2 - 2s^k^ - 2k^S^ + k^) 


The derivafives of Y(k so is (t>xo yo) respecf fo fhe confrol paramefers can be found by 
sfraighfforward calculafions using fhe above. 
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